1 Purpose

This document summarizes preliminary results of analyses that seek to quantify the extent to which salmon survival (log(R/S)) is correlated with (a) ocean conditions during early marine life and (b) potential salmon competitors later in marine life, and the degree to which these relationships vary over time and space. More details on motivation for project can be found here and all code to reproduce the analyses in this document is here. All findings are preliminary and subject to change.

2 Methods

Our overarching approach is to fit four general classes of spawner-recruitment models to characterize relationships between ocean conditions and salmon productivity over space and time:

  1. Stationary models that estimate invariant relationships between ocean conditions and salmon productivity (e.g., Connors et al. 2020)

  2. A-priori change point (or “Era”) models that allow relationships to vary among periods that corresponded with observed Northeast Pacific regime shifts in 1976/77 and 1988/89 (e.g., Malick 2020), as well as the onset of the recent marine heat wave period in 2014 (in upcoming analyses)

  3. Random walk models that allow relationships to evolve gradually through time (e.g., Malick 2020)

  4. Hidden Markov models that allow relationships to vary according to the state (regime) the system is in, where both the sequence of states and state-specific relationships are estimated using a Hidden Markov Model framework.

Details on model formulation are provided below. These models are fit in a Bayesian estimation framework uand inference is primarily based on the magnitude, direction, and uncertainty of standardized covariate effects and the degree to which they vary over time and across ocean regions, species, and life histories.

2.1 Data

2.1.1 Spawner-recruitment

We have been working to update and add to a large collection of time series of spawner (escapement) and total recruitment (catch plus escapement) for all five species of Pacific salmon. Details on the specific timeseries are available here. The spatial and temporal extent of those populations and species that we have fit models to so far are detailed below.

Figure 1. Survival (log[R/S]) time series of all Sockeye stocks included in analyses. Vertical dashed lines correspond with pre-defined ocean regime shifts included in the “eras” model.

2.1.2 Covariates

Our analysis is focused on two covariates:

  1. Sea surface temperatures near juvenile salmon ocean-entry points to index regional-scale environmental variability during the first few months in the ocean, which are hypothesized to represent a critical period in the Pacific salmon life cycle that can strongly influence year class strength (Mueter et al. 2002).

  2. The abundance of pink salmon across the North Pacific in the second and/or third years of marine life as an index of potential direct or indirect competition for food. This approach is consistent with research that has suggested some salmon species primarily exhibit responses to competition with other salmon during their second and/or third growing seasons at sea (Connors et al. 2020; Ruggerone et al.2023).

Figure 2. Sockeye ocean entry locations and associated climate and competitor time series. Numbered stocks are named in Figure 1. (a) Unique ocean entry points of Sockeye stocks in four Ocean Regions defined by ocean entry points. (b) Timeseries of Pink salmon abundance index in millions (top row) and SST anomaly (bottom row). Timeseries are age-weighted to reflect average ocean conditions experienced by brood cohorts regardless of age at ocean entry. Stock-specific timeseries are shown by light lines and regional averages by bold lines.

2.2 Stationary models

In the simplest stationary class of models, salmon productivity was modelled hierarchically as a function of spawner abundance, early marine ocean conditions, and potential competitor abundance later in marine life:

\[\begin{equation} ln(R_{i,t}/S_{i,t}) = \alpha_{i} - \beta S_{i,t} + \gamma_{k,i,r} X_{k,i,t} + \varepsilon_{i,t} \tag{2.1} \end{equation}\]

where \(R\) is the total recruitment from population \(i\) and spawners \(S\) in brood year \(t\), \(\alpha\) is productivity (intrinsic rate of growth), \(\beta\) is the magnitude of within brood-year density-dependent effects on survival, \(\gamma\) is the effect of covariate \(k\), \(X\) is the specific covariate standardized to mean 0 and standard deviation of 1, and \(\epsilon\) is inter-annual variation in survival which was assumed to follow a first-order autocorrelated process:

\[\begin{equation} \varepsilon_{i,t}=\phi \varepsilon_{i,t-1} + \sqrt{1-\phi^2 }\delta_{i,t} \\ \delta_{i,t}\sim N(0,\sigma_{i}^2) \tag{2.2} \end{equation}\]

where \(\phi\) is the autocorrelation coefficient, and \(\delta\) represents uncorrelated, white noise, interannual variation in survival.

The population-specific parameters \(\gamma_{k,i,r}\) were modelled hierarchically by assuming they arise from common hyper-distributions:

\[\begin{equation} \gamma_{k,i,r}\sim N(\mu_{\gamma,k,r},\sigma_{\gamma,k,r}^2) \tag{2.3} \end{equation}\]

where \(r\) is one of four broad ocean regions within which salmon populations have previously been shown to exhibit similar responses to ocean climate and ecosystem variation and \(\mu\) is the region level average effect.

For Fraser River Sockeye stocks (No. 1-19 in Fig. 1), spawner abundances were in units of effective female spawners (i.e., female spawner abundance adjusted for unspawned eggs), whereas for all other stocks spawner abundances were total male and female spawners. Thus, \(\alpha_{i}\) parameters were split into two groups (one group that included Fraser River stocks and another group that included all non-Fraser River stocks) that were exchangeable within each group but not between the two groups.

2.3 Era models

The first class of non-stationary models we considered allowed the \(\gamma_{k,i,r,e}\) parameters to change at predefined points in time corresponding with well documented Northeast Pacific regime shifts in 1976/77 and 1989/99 (Hare and Mantua, 2000) resulting in three eras (\(e\)). As with the stationary models, we assumed the era-specific climate and ecosystem parameters were exchangeable among populations within the same ocean ecosystem and era, thereby estimating both population specific and ocean ecosystem level effects by era.

2.4 Random walk models

The second class of non-stationary models we considered allowed the climate and ecosystem parameters to evolve over time according to a random walk:

\[\begin{equation} \gamma_{k,i,t}=\gamma_{k,i,t-1}+w_{k,t}\\ w_{k,t}\sim N(0,\sigma_{w_{k}}^2) \tag{2.4} \end{equation}\]

where \(w\) is process variation and \(\sigma_{w}^2\) determines the degree of temporal variability in the \(\gamma\) series. For these models the covariate effects were modelled independently instead of being exchangeable among populations within regions (i.e., hierarchically).

2.5 Hidden Markov models

The last class of non-stationary models is conceptually similar to the random walk, but where the coefficient effects are pooled into discrete ordered distributions (i.e. low or high \(\gamma\)) representing unobserved regime states \(q\), and the probability of being in each state is evaluated at each time-step. The estimated timeseries of regime states \(z_{t}\) is modelled as a Markov process (ie. year-to-year regime states depend only on the state in the previous time-step (\(p(z_{t}|z_{t-1})\))) that progresses according to an estimated transition probability matrix.

Covariate effects \(\gamma_{k,i,q}\) are modelled as independent among stocks, not pooled by ocean ecosystem. The model requires that one primary covariate (SST) defines regime states, and the other (Competitors) has discrete \(\gamma\) distributions corresponding with each state.

2.6 Model fitting

We fit the model described above in a Bayesian estimation framework with Stan Stan Development Team 2020, which implements the No-U-Turn Hamiltonian Markov chain Monte Carlo algorithm Hoffman et al. 2014 for Bayesian statistical inference to generate the joint posterior probability distribution of all unknowns in the model. For each model run, we sampled from 4 chains with 4,000 iterations each and discarded the first 1,000 as warm-up. We assessed chain convergence visually via trace plots and posterior predictive checks and by ensuring that \(\hat{R}\) [potential scale reduction factor; Vehtari et al. 2021] was less than 1.01 and that the effective sample size was greater than 400 for each parameter in the model.

3 Results

3.1 Preliminary Insights

3.1.1 Latitudinal responses

We found strong evidence for latitudinal responses to both ocean warming and interspecific competition, revealed by the stationary models (Figure 3). The effect of SST ranges from moderately negative in the southernmost region to strongly positive in the northernmost region. Conversely, competitor abundance has a strongly negative effect in the south but a weak effect in the north (Figure 3).

3.1.2 Evidence for nonstationarity

The various classes of nonstationary models suggest there is evidence for variably strengthening, weakening, or emerging relationships over time and across regions. For example, assuming predefined ocean regimes in the Era model resulted in a negative, slightly shifting SST relationship in the West Coast region (Figure 5), while the Random Walk model suggests a gradually declining relationship from slightly positive in the earliest years (Figure 7). Across Alaska regions, the era model predicts that the 3 northern regions similarly shift from weakly positive to slightly negative relationships with competitors (Figure 5), but the random walk model predicts relatively stationary relationships over time in Southeast Alaska and the Gulf of Alaska (Figure 7).

The Southeast Alaska region may exhibit unique relationships due to its ocean conditions being highly dependent on year-to-year shifts in the North Pacific current. While the era model estimates relatively weak SST effects in Southeast Alaska across time periods (Figure 5), the random walk model suggests directional change that is dramatically negative in the last two decades (Figure 7).

Some common insights arise from the era and random walk models, including:

  • A strong negative shift in response to competitor abundance in the Bering Sea, where Pink salmon abundance more than doubles over the time period examined (Figure 2);

  • A strengthening of the already negative relationship of West Coast stocks with competitor abundance;

  • The most overall stable dynamics occurring in the Gulf of Alaska.

The Hidden Markov models show relatively weak evidence for non-stationarity. In general, the HMM predicted that covariate relationships readily alternated between states, rather than exhibiting periods of stationarity separated by regime shifts like we would expect this type of model to detect. As an exception, the bi-directional pattern of the SST effect in Southeast Alaska is also visible in the HMM results (Figure 8).

3.2 Stationary models

Figure 3. Posterior probability distributions of the predicted effect of SST (top), competitors (middle), and the combined effect (bottom) from all covariate terms, on Sockeye salmon survival. Regional hyper-distributions of the covariate effects are in bold lines, with individual stock-specific distributions illustrated by the light lines. Covariate effects are standardized (i.e., per standard deviation unit increase in each covariate).

Figure 4. Stock-specific posterior mean (points) and 95% CI (horizontal bars) estimates of SST and Competitor covariate effects (left and right, respectively). Vertical lines and shaded areas are regional hyperdistribution means and 95% CI.

3.3 Era models

Figure 5. Posterior probability distributions of the predicted effect of SST and competitors on Sockeye survival over three pre-defined time periods/eras (earliest in top panel). Breakpoints between eras correspond with documented North Pacific ocean regime shifts. Regional mean (hyper-distribution) covariate effects are in bold lines and individual stocks’ distributions are in light lines.

Figure 6. Stock-specific posterior mean (points) estimates of SST and Competitor covariate effects through three pre-defined eras: <1976 (lightest shades), 1977-1988 (middle shades), and >1989 (darkest shades). Vertical lines and shaded areas are regional hyperdistribution means and 80% CI.

3.4 Random walk models

Figure 7. Time-varying posterior mean estimates of SST and Competitor covariate effects, modelled as a random walk. Individual stock estimates are in faint lines, while regional means and 90% CI are represented by bold lines and shaded areas. Regional summaries are post-hoc calculations, as random walk \(\gamma\) parameters are not grouped hierarchically.

3.5 Hidden Markov models

Figure 8. Time-varying mean estimates of SST and Competitor covariate effects from hidden Markov models. Time-varying effects are calculated from the probability timeseries of regime states \(p(z_{t})\) and covariate coefficients for each state \(\gamma_{k,i,q}\) (Figure 9). Individual stock estimates are in faint lines, while regional means and 90% CI are represented by bold lines and shaded areas. Regional summaries are post-hoc calculations, as each stock was modelled separately.

Figure 9. Posterior estimates of SST and Competitor coefficients \(\gamma_{k,i,q}\) for each hidden state \(q\) in hidden Markov models. Vertical ticks show posterior means, shaded areas are 80% CI, and horizontal bars are 95% CI. States are determined by SST, and a competitor effect is then determined for each state.

4 Next steps

4.1 Data and environmental covariates

  • Update Sockeye, Pink, and Chum spawner-recruitment time-series based on data from ADF&G and DFO; incorporate into analysis (in progress)

  • Incorporate Chinook and coho spawner-recruitment time series into analysis

  • Explore alternative temporal and spatial domains over which to average SST anomalies to derive the early ocean temperature covariate

  • Explore alternative indices of interspecific competition at sea, including combined abundance of pink, chum, and Sockeye, as well as effective density to account for body size

4.2 Modelling

  • Explore hierarchical parameterization of random walk model (and hidden markov?) to estimate region-wide effects

  • Add a ‘marine heatwave’ time period to era models that begins in ~2014 at the onset of several anomalous warm years at sea